function [log_likelihood,nan_pct]=ml_time_score_types_cdf(theta)

%%% calculates the log-likelihood allowing multiple factors. Takes
%%% simulated factors and risk aversion from step1j as input 



 global X_ind_s scores ind_id factor_n  ind_n draw_n  time_choice_n time_choice_ind  n_type sim_theta_i sim_kappa_i est_p_type  sim_1j_contribs maxima minima  r_scale sig_th_lower_limit CRRA2  K_time est_kappa_vector CARA x1 t1 x2 t2 hyper analytical fine_grid
    



    
    %%% simulated factors from step1j
    sim_factors=scores;
    
    %%% Probability contribution from previous step (likelihood of measures
    %%% and lottery choices
    %%% observed for each individual given his characteristics and
    %%% simulated draws)
    sim_l_contrib_1j=sim_1j_contribs;
    
    % # of observations including rows with simulated factor draws for each
    % person
    rows=size(sim_factors,1);
    
    sex=X_ind_s(:,1);

    %%% simulated risk aversion from step1j given characteristics and
    %%% simulated factors  
    theta_i=sim_theta_i;
    
    %%% taken from previous iteration where kappa is estimated in the risk
    %%% phase 
    if K_time==0
        kappa_i=sim_kappa_i;
    end
    
    %%% For robustness with a separate mistake intercept for intertemporal
    %%% tasks
    if K_time==1
        kappa_types=est_kappa_vector(1:n_type);
        kappa_coeffs=est_kappa_vector(n_type+1:end);
        %%% adds the extra time intercept for kappa to the end of the
        %%% coefficient list
        kappa_extra=theta(end);
        theta(end)=[];
    end    
    
    %%% estimated type prevalences from step1j
    p_type=est_p_type;
   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Time Part %%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    %%% parameter assignment for time analysis
    
    r1=theta(1);
    r2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[];        

    sig_r1=theta(1);
    sig_r2_n=theta(2:factor_n+1);
    theta(1:factor_n+1)=[]; 
    
    type4_5=NaN(2,n_type);
    
    sim_l_contrib_time=NaN(ind_n,time_choice_n);
    av_l_part=NaN(ind_n,n_type);  

    %%% weigh over unobserved types, each has its own set of intercepts for
    %%% the structural parameters of interest
    for tt=1:n_type
        %%% last type's proba is just 1-sum(other probas) so don't calculate
        %%% type intercepts
        type4_5(:,tt)=theta(1:2);
        theta(1:2)=[]; 
        
        if K_time==1
            kappa_i(:,tt)=0.5*normcdf(kappa_extra+kappa_types(tt)+sex*kappa_coeffs(1)+sim_factors*kappa_coeffs(2:end));
            kappa_i=min(maxima(2), kappa_i);
            kappa_i=max(minima(2), kappa_i);
        end
        
        %%% individual delay aversion based on sex and particular factor draws -
        %%% deterministic part - note that as above, the factor is determined
        %%% by the orthogonal draw AND by X*alpha - a function of observables

        %%% Even though the following parameters are individual-specific - thus
        %%% _i - they vary here within individuals due to the different factor
        %%% draws used in each row for the factor simulation
        r_i=r_scale*normcdf(type4_5(1,tt)+sex*r1+sim_factors*r2_n);

        %%% individual variability of delay aversion based on sex and particular
        %%% factor draws
        sig_r_i=normcdf(type4_5(2,tt)+sex*sig_r1+sim_factors*sig_r2_n)+sig_th_lower_limit;

        
        r_i=min(maxima(4), r_i);
        r_i=max(minima(4), r_i);
        sig_r_i=min(maxima(5), sig_r_i);
        sig_r_i=max(minima(5), sig_r_i);       

        if analytical==1 
                %%% analytical calculation of indifference thresholds for
                %%% intertemporal choice tasks

                %%% first what utility to consider (CRRA, CRRA2, CARA)
                if CRRA2==1
                    %%% No longer need to limit risk aversion at 1 as with standard CRRA to keep
                    %%% utility positive for discount rate calculations. At
                    %%% 10, indifference thresholds are zero  
                    theta_limit=min(theta_i(:,tt), 10);   
                    theta_limit=max(theta_limit, -.3);            
                    
                    %%% utility of either option in period in which it is
                    %%% received
                    U_1=((x1.^(1-theta_limit)-1)./(1-theta_limit));  
                    U_2=((x2.^(1-theta_limit)-1)./(1-theta_limit));

                    %%% as with standrd CRRA, at theta=1, U=log(X1)
                    [xxx,yyy]=find(theta_limit==1);
                    U_1(xxx,:)=repmat(log(x1),size(xxx,1),1);  
                    U_2(xxx,:)=repmat(log(x2),size(xxx,1),1);   
                elseif CARA==1
                    %%% as above, relaxes upper limit on theta for discount
                    %%% rate calculation as CARA stays 
                    %%% positive for positive $ amounts.
                    theta_limit=min(theta_i(:,tt), 1);   
                    theta_limit=max(theta_limit, -.01);            

                    U_1=(1-exp(-theta_limit.*x1))./theta_limit;   
                    U_2=(1-exp(-theta_limit.*x2))./theta_limit;

                    %%% for CARA, at theta=0, U=x1
                    [xxx,yyy]=find(theta_limit==0);
                    U_1(xxx,:)=repmat(x1,size(xxx,1),1);  
                    U_2(xxx,:)=repmat(x2,size(xxx,1),1);               
                else
                    %%% Need to keep utility positive, so theta<1
                    theta_limit=min(theta_i(:,tt), .99);
                    theta_limit=max(theta_limit, -.3);            
                    U_1=((x1.^(1-theta_limit))./(1-theta_limit));
                    U_2=((x2.^(1-theta_limit))./(1-theta_limit));
                end

                %%% second, what discounting to consider: exponential/hyperbolic
                if hyper==1
                    %%% indifference threshold given task parameters and the individual's
                    %%% estimated coefficient of risk aversion
                    r_xy_comp=((U_2./U_1)-1)./(t2-t1);    
                else
                   %%% indifference threshold given task parameters and the individual's
                   %%% estimated coefficient of risk aversion                  
                    r_xy_comp=(U_2./U_1).^(1./(t2-t1))-1;
                end

        end



        %%% Assumes lognormal errors on the delay aversion parameter, so log:
        %%% ln(r/sig_r)~N(0,1).      
        r_i2=r_i.^2;
        %%% variance of the error term (with a lognormal distrib) on the discount
        %%% rate
        var=(sig_r_i).^2;
        %%% Formulas for mean and sd of normal distrib of which the lognormal
        %%% distrib with mean=r_i and sd=sig_r is the exponent. Done to get pdf,
        %%% cdf of the lognormal distrib which in Matlab takes as parameters the
        %%% mean and sd of the corresponding normal distrib of which it is the log.
        % mean of standard normal distrib
        mu=log(r_i2./sqrt(var+r_i2));
        % sd of standard normal distrib. Adds the small number at the end
        % for purposes of numerical optimization
        sigma=sqrt(log(var./r_i2+1))+(1*10^-20);

        mu_full_time=repmat(mu,1,time_choice_n);
        sigma_full_time=repmat(sigma,1,time_choice_n);

        %%% Probability of choosing the later option i.e. probability that the
        %%% individual's discount rate is below the threshold for a
        %%% particular choice task, given the individual's coefficient of
        %%% risk aversion. Need to take the log here of the 
        %%% resulting difference divided by sig_r as it is the log of the
        %%% error which is normally distributed. 
        p_distant=logncdf(r_xy_comp,mu_full_time,sigma_full_time);
        p_near=1-p_distant;

        % kappa percent of the time and individual makes the wrong choice.
        choice_distant=zeros(size(p_distant));
        for i=1:time_choice_n
            choice_distant(:,i)=p_distant(:,i).*(1-kappa_i(:,tt))+p_near(:,i).*kappa_i(:,tt);
        end
        choice_near=1-choice_distant;
        sim_l_contrib_time=(choice_distant.^time_choice_ind).*(choice_near.^(1-time_choice_ind));   
   
    
    
        
        %%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions from the various
        %%%%%%%%%%%%%%%%%%%%%% types of measures and lottery and
        %%%%%%%%%%%%%%%%%%%%%% intertemporal task choices 
       
        % for numerical estimation purposes
        sim_l_contrib_time_plus=sim_l_contrib_time+1*10^(-20);

        sim_l_contrib_total=[sim_l_contrib_1j(:,tt) sim_l_contrib_time_plus];

        % Individual likelihood contribution conditional on factor draws
        sim_ind_draw_contrib=prod(sim_l_contrib_total,2);

        % individual likelihood of measures and choices averaged
        % over simulated factor draws
        av_l_contrib=cell2mat(accumarray(ind_id,sim_ind_draw_contrib,[],@(x){sum(x)}))/draw_n;
    
        % In order to avoid zeros for the log
        av_l_part(:,tt)=av_l_contrib+1*10^(-200);
        
        sim_l_contrib_time=NaN(ind_n,time_choice_n);
        clear sim_l_contrib_total sim_l_contrib_total_plus sim_ind_draw_contrib av_l_contrib r_i sig_r_i theta_i_cat r_xy_comp r_i2 var mu sigma mu_full_time sigma_full_time p_distant p_near choice_distant choice_near

    end

    
    %%%%%%%%%%%%%%%%%%%%%%% Combining the likelihood contributions of the
    %%%%%%%%%%%%%%%%%%%%%%% various unobserved types and weigthing them by each type's
    %%%%%%%%%%%%%%%%%%%%%%% estimated population prevalence
    
    av_l=zeros(ind_n,1);
    for tt=1:n_type
        av_l=av_l+p_type(tt)*av_l_part(:,tt);
    end
    
    log_likelihood=-nanmean(log(av_l));
    
    %%% will capture the % of nan and inf in the final average likelihoods
    nan_pct=1-mean(isfinite(av_l));
    if nan_pct>.05
        log_likelihood=NaN;
    end
    
end